Symmetry restoration in Hartree-Fock-Bogoliubov based theories 
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We present a pfaffian formula for projection and symmetry restoration for wave functions of the 
general Bogoliubov form, including quasiparticle excited states and linear combinations of them. 
This solves a long-standing problem in calculating states of good symmetry, arising from the sign 
ambiguity of the commonly used determinant formula. A simple example is given of projecting good 
particle number and angular momentum from a Bogoliubov wave function in the Fock space of a 
single j-shell. 



Introduction. The Bogoliubov transformation offers a powerful way to introduce correlations into multi-fermion 
wave functions. The variational theory based on it, the Hartree-Fock-Bogoliubov (HFB) theory, has been very useful 
^"*^ ■ in nuclear physics. However, the variational wave functions need not respect symmetries of the Hamiltonian, hindering 
its use for spectroscopic purposes. An obvious fix is to project the wave functions onto eigenstates of the conserved 
quantum numbers. However, present methods to carry out the projection are beset with technical difficulties. The 
purpose of this letter is to present a projection formula that is applicable to general Bogoliubov wave functions, 
including those for odd particle number. The results are generalized for the evaluation of overlaps, as those required 
in configuration mixing theories based on HFB wave functions, commonly referred to as Generator Coordinate Method 
(GCM). 

We first remind the reader that an operator Vk for projecting onto a symmetry group representation K is given 
by the integral 



VKi = ^JdnRf i (si)K(n), (i) 



(N 
> 

' Here R is an operator of the symmetry group, is a diagonal element of a matrix representation of the group, dx 
is the dimension of the representation matrix R K and tto = J dil is the volume integral over the group. The main 
conserved quantum numbers that we wish to restore in nuclear physics are particle number N and angular momentum 
J. These are both very familiar but for concreteness we note that particle number is associated with the gauge group 
U(l) and the group integral is d<j> where is the gauge angle. In the case of angular momentum, the integration 
is over the Euler angles sin(/3)e?a(i/3d7 and the representation matrices are the Wigner D-functions. The probability 
of the component with quantum number K in the state \w) is given by the integral 

(K\w) 2 = (w\T K \w) = ^ J <KIR$ (w\K\w). (2) 



In this letter we treat only the problem of calculating the overlaps; for applications one also needs to calculate 
matrix elements of physical operators. In the past, the computation of the overlap (w\lZ\w) was carried out by the 
Onishi formula[lJ (see also 0, Eq. E.49]). Unfortunately, the formula has square root sign ambiguity which makes 
it useless for projection, except in special cases. Several suggestions have been made in the past for overcoming this 
sign problem 3|-|6[ . In ref. [6] , Robledo proposed a promising new formula based on the pfaffian rather than the 
determinant. However, his formula requires the inverse of the Bogoliubov transformation matrix U, and is thus not 
applicable to wave functions for which the U matrix is singular. This is the case for all wave functions that have zero 
overlap with the vacuum. In particular, the formula cannot be used directly for states of odd particle number. 

Here we propose a pfaffian expression which can be easily extended to odd-iV wave functions, and indeed to states 
with more than one quasiparticle excitation. To establish the notation, we write the effect of the symmetry operation 
as 

TZcjlZ- 1 = Rijc] ; TlciRr 1 = R *j c J ( 3 ) 

3 3 

where cj and Cj are the usual Fock-space creation and annihilation operator in some convenient basis. Note that the 
matrix R depends on the specific details of the basis states and does not have to belong to an irreducible representation 
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of the group. The wave function is characterized by the U, V matrices of the Bogoliubov transformation. Use is made 
of the Bloch- Messiah decomposition (see [2[ for details and notation) that expresses those matrices as the product of 
unitary matrices D and C and special block diagonal matrices U and V, namely U = DVC and V = D*VC @, Eq. 
7.8]. The unitary D transformation defines the "canonical" basis with creation and annihilation operators and a. 

We first consider the simpler case in which the wave function has a non-zero overlap with the vacuum. Then it can 
be expressed in the canonical basis as 



Y[(u a + v a alal)\). 



(4) 



Here n is the number of pairs in the wave function and the matrices U, V have dimension (2n x 2n). To specify the 
phase of the wave function, we may take all u a positive definite. The overlap in this case is given by 



(w\K\w) 



(-1)" 

n>«) 
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V T U V T R T V* 
-V^RV WV* 



(5) 



where pf (M) is the pfaffian of the matrix M. We outline an alternative derivation below. Note that to use Eq. §5§ 
the U, V matrices in the canonical basis must be truncated to omit columns for which v a — (see also Ref. [H). This 
simply means omitting the part of the Fock space that is not occupied. 

The generalization of Eq. ([5]) to deal with arbitrary overlaps is straightforward but it requires to write the wave 
function \w) of Eq. ((4|) as 



detC 

H = TT" — PlP2 • 



n 



..&»i> 



(6) 



a=l u a 



where the /3 M are Bogoliubov quasiparticle annihilation operators with amplitudes U and V. Here an unnormalized 
wave function is obtained by the product of all the Bogoliubov-transformed annihilation operators acting on the 
vacuum, and [] , v~ is the normalization factor. The phase det C is required for consistency of all the definitions. 
The overlap is then given by 



, , . , det C* det C" 

(w\K\w ) = (-1) ~ pf 



V T U V T R T V 
-V'^RV WV* 



(7) 



This formula is useful in dealing with configuration mixing of symmetry restored HFB wave functions, as required in 
the implementation of the most general version of the Generator Coordinate Method (GCM). The connection between 
Eq. ([7]) and Eq. (7) of @ is not straightforward and requires of some lengthy calculations, details are given in 0]. 

Eq. (JTJ) may be extended to wave functions that are orthogonal to the vacuum by considering the more general 
canonical form 



(8) 



Again, if the canonical U, V (U' , V) matrices are truncated to omit columns for which v a = (v' a = 0), Eq. (|7|) is 
still applicable. In the case of odd- TV ground states, only a single additional operator is needed, 



\qw 



(9) 



We use the notation q for the row vector of the coefficients qi (qi^ in matrix notation), and for the row vector of 
zeros. Then the generalization of Eq. is 



{qw\Kqw = -1 — = — pf 

1L v a < 



V T U T 


-q'RV -q'i?qt 

-V'^RV -F'tflqt 



V T R T c{ T V T R T V* 
q*R T q' T q*R T V* 





T 







(10) 



The shape of this matrix is (2n + 2 + 2n) x (2n + 2 + 2n). To derive it and Eq. ([7]), first note that the expectation 
value of a product of single- fermion operators on is given by the pfaffian of all possible contractions (il. Hoj 



(axa 2 . ..a k ) = pf ( Si j ) 



(11) 
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where Sij is a skew-symmetric matrix with upper triangular elements Sij — (\aiOtj\) (i < j) The overlaps can be 
written in this form using Eq. ([6]). The overlap in Eq. (| 10[) is derived by evaluating the contractions in the operator 
product 

(PlAn-X---Pl^ q ~P'i~&---K) (12) 

where j3' = Hfi'llr 1 etc. The matrices entering the pfaffian of Eq. (fT0|) are easily identified: 

(V T U),„ = (|4/3t|) (13) 

(v T R T q ' T ), = (144,1) = J2^M^\) ( 14 ) 

3 

(V T R T V'*)„ = (1441) (15) 
q*i? T q' T = J^q^^c],]) (16) 

33' 

(q*R T V>*)„ = £<Z*(M:i> (17) 

3 

(crvv = (i44i). (is) 

The generalization of Eq. (|10[) to multi-quasiparticle overlaps, with r annihilation operators /3 Mj . (Bogoliubov 
amplitudes U, V) to the left of TZ and s creation operators j3' v . to the right, is tedious but straightforward 



(w% r ... •••4tK> - (-l)"(-l)Kr-D/2 detC*det^ pf 



a 



U T t/ F T pt T/ T i? T q' T V T R T V* 

p*V q*pt q*i? T q' T q*i? T y'* 

-q'-RF -q'i?qt p'q' T pT" 

-y'ti?y -yt^qt -y'V 7 t/'V* 



(19) 

For this expression to make sense both r and s must have the same number parity. The objects p and q (p' and 
q') are matrices of dimension r x 2n (s x 2n) with matrix elements p Mjm = V mfJlj and q^m = U mtij . If some of 
the /3 annihilation operators are replaced by creation ones fft the appropriate rows of q and q have to be redefined 
accordingly. It is easy to check that Eq. (|T5]) reduces to Eq. (TTU)) in the limit p = 0. Apart from the fact that Eq. 
(I19p includes the phase of the matrix element, this expression has the advantage over more traditional approaches 
[l2| that the combinatorial explosion in the evaluation of the left hand side of Eq. (fill)) , namely the fact that (r + s)!! 
contractions have to be considered if the multi-quasiparticle overlap is computed with the standard Generalized Wick's 
theorem, is completely avoided (see |11| for another approach based on the finite temperature formalism). 

Exam/pie. As a proof of principle, we carry out the projection for an odd-iV wave function having a nontrivial 
structure with respect to angular momentum and particle number. We take the Fock space as the 6-dimcnsional 
space of orbitals in a j — 5/2 shell. The creation operators are labeled by azimuthal angular momentum j z = m. 
The wave function for the test is 

\qw) = c\ /2 (u + «4 /2 cl 5/2 ) |) (20) 

with (u, v) = (0.8,0.6). We project simultaneously on particle number and angular momentum with the operator 
'Pn'Pj.Jz ■ We use a 4-point uniform mesh for integrating the gauge angle and a 5-point Gauss-Legendre quadrature for 
integrating over the angular variable cos(/3). There is no necessity to integrate over the other Euler angles because the 
wave function Eq. ([20]) is an eigenstate of J z . The results are shown in Table I. The projected quantum numbers N 
and J are given in the first two columns. The third column gives the exact decomposition, and the fourth column the 
numerical projection. One sees that there is complete agreement to the level of numerical precision in the integrations. 

Discussion. Besides the overlap function (w\lZ\w) we need the matrix elements of various operators in the symmetry- 
restored states. For the most important operators they can be expressed as a single integral over matrix elements 
of the type (w\01Z\w) where O is an operator such as the Hamiltonian. It is straightforward to calculate this 



operator matrix elements using the Balian-Brezin formula[12j or the multi-quasiparticle overlap of Eq. (|19|) . Unlike 



the formulation in Ref. our method here does not require one to construct quasiparticle states explicitly. Our 
procedure is easily extended to multi-quasiparticle matrix elements with a final result that avoids the combinatorial 
explosion that plagues other methods used to evaluate those overlaps. For an fc-quasiparticle excitation, the pfaffian 
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TABLE I: Test of number and angular momentum projection for the wave function of Eq. (|20[) . for which J z = 1/2. 





(NJJ z \qw) 2 


N J 


analytic 


numerical 


1 3/2 





0.00000 


1 5/2 


u 2 = 0.64 


0.64000 


3 1/2 





0.00000 


3 3/2 


v 2 /7 « 0.05143 


0.05143 


3 5/2 


u 2 /2 = 0.18 


0.18000 


3 7/2 





0.00000 


3 9/2 


5« 2 /14 ~ 0.12857 


0.12857 



matrix is augmented by 2k rows and columns. A program that demonstrates the method for the example in Table I 
is provided in the supplementary material. After this work was posted on arXiv 11311 we learned that a similar fully 
general formula for the overlap was obtained independently by Avez and Bender |l4| and Oi and Mizusaki [15| . 
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Supplementary material 

We want to prove the equivalence between the result of Eq. (7) in the letter for the special case TZ = 1 

, , ,s , 1 Nndet(C*)det(C) , / V T U V T V* \ 

(w\w) = (-1)™ Wrr r — -P f ,+ , + , • ( A - 21 ) 

V 1 1 V ' \-V'W U«V* J 

and the overlap formula of Eq (7) in Ref. @. In both cases U, V are the 2n x 2n Bogoliubov transformation matrices 
and C the third (unitary) transformations of their Bloch-Messiah decomposition. The normalized wave function \w) 
is written as 

H = ^l») (A.22) 
with \w) = Pi . . . /02n|}- Using the Bloch-Messiah decomposition it is easy to show that \w) and the wave function 



= exp ( 2^ Mkk '°* C V 

\ kk' , 



of Eq (1) in |6| are related by 

|0}=pf(CTV)W (A.23) 

In terms of \w) 

/ -,\nvf vTu V T V'*\ ... 

lw\w'} = (-l) n pf ,, , . (A.24) 

v \ -v'^v u'W* ) 

In Ref. [f|, Eq (7) the following equation is derived for the overlap. 
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where we write M = (VU^ 1 )* and N in Eq (7) of 6] is twice the n here (N = 2n, and therefore S N = (-l) N ( N + 1 )/ 2 = 
(— 1)"). To prove the equivalence of the two expressions, we make multiple use of the pfafhan identity for the 
congruence transformation, 

tf(A T BA) = det(A) pf(B). (A.26) 

We first write the matrix in Eq. (|A.24[) as 



V T U V T V>* \ ( V T UV~ X I \ V 
-yty [7'ty* J _ I o V'i I \ -I -u^V*- 1 j \ V* 



(A.27) 



Then using Eq. (|A.26[) the overlap in Eq. (|A.24[) can be expressed as 



<*K =(-i)"«i<M(r)d«M ( rMpi ( M V _ J ^ / _ 1 ) . (A.28) 



We next use the congruence transformation to express the matrices in Eq. (|A.28I) and (IA.25|) as 
( M*- 1 I \_( 1 \ f M*- 1 \ (l M*\ 

y -i -M'- 1 ) ~ \ -m* I J I -Af'- 1 + m* y y I i 

/ M' -I \ _ / I o\ f M' \ / I —M'~ l \ 

\ I -A/* / _ 1 Af'- 1 I J I -M* + M'~ l ) \ I / ' 



(A.29) 



(A.30) 



The congruence transformation matrices have determinant equal to one, so the pfaffians of left-hand side are are equal 
to those of the transformed matrices, pf(M* - 1 )pf(-M /_1 + M*) and pf (M')pf (M' -1 — M*), respectively. With these 
results and the properties pf(A _1 ) = (-l)"/pf(A) and pf(-A) = (-l) n pf(yl) we get 

and 

= pf(A/')pf(-Af'- 1 + AT) (A.32) 
The equivalence between the two expressions Eqs. (|A.31[) and (|A.32|) is evidenced by taking into account that 

^1 = det(V T ) P f(V T -iU T ) = P f(^F) (A.33) 

which finishes the proof. 

The generalization to the case where 72. is not the unity matrix, proceeds along the same lines. 
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